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Abstract 

A new force balance model for the EFIT [l[ magnetohydrodynamic equilibrium technique for tokamaks is 
presented which includes the full toroidal flow and anisotropy changes to the Grad-Shafranov equation. The 
free functions are one dimensional poloidal flux functions and all non-linear contributions to the toroidal 
current density are treated iteratively. The parallel heat flow approximation chosen for the model is that 
parallel temperature is a flux function and that both parallel and perpendicular pressures may be described 
using parallel and perpendicular temperatures. This choice for the fluid thermodynamics has been shown 
elsewhere 0] to be the same as a guiding centre kinetic solution of the same problem under the same 
assumptions. The model reduces identically to the static and isotropic Grad-Shafranov equation in the 
appropriate limit as different flux functions are set to zero. An analytical solution based on a modified 
Soloviev solution for non-zero toroidal flow and anisotropy is also presented. 

The force balance model has been demonstrated in the code EFIT TENSOR by modifying an existing 
code EFIT++ Q. Benchmark results for EFIT TENSOR are presented and the more complicated force 
balance model is found to converge to force balance similarly to the usual EFIT model and with comparable 
speed. 

Keywords: tokamak, equilibrium, EFIT, reconstruction, flow, anisotropy, MHD, magnetohydrodynamics, 
nuclear, fusion, plasma, fast particles, NBI, ICRH 



1. Introduction 

The macroscopic equations of magnetohydrody- 
namics (MHD) provide the basic starting point for 
an understanding of plasma physics in a modern 
tokamak experiment. Good knowledge of the to- 
tal equilibrium force balance provides information 
about the magnetic topology and the plasma ther- 
modynamic variables so that more detailed stability 
and transport treatments can be pursued. The suc- 
cess of the tokamak concept has been possible, in 
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part, due to the simple models which support infer- 
ence of the plasma configuration from incomplete 
measurements of related parameters. Most basic 
physical variables can only be measured indirectly 
on a fusion experiment through a sophisticated and 
expensive set of diagnostics. Good physical models 
of the configuration are essential, particularly for 
future power stations which cannot accommodate 
complex diagnostics programmes. 

Modern tokamak experiments contain a signifi- 
cant portion of fast ions [J] resulting from heat- 
ing processes such as neutral beam injection (NBI) 
and ion-cyclotron resonance heating (ICRH) which 
can rotate the plasma and also produce highly 
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anisotropic fast particle pressures 0, [fjj. Both of 
these effects can move the plasma off magnetic 
field surfaces and significantly alter the density 
profile and magnetic topology (see, for example, 

The tokamak equilibrium reconstruction code 
EFIT [lj has served as the de-facto standard tech- 
nique to infer equilibrium from experimental diag- 
nostics and there have been many different code 
implementations of this technique. EFIT solves the 
MHD force balance for a static and isotropic pres- 
sure, although approximate inclusions of toroidal 
flow Q and anisotropy 11 1 exist in some versions 
through modifications to pressure function. 

In this paper, we describe an extension of the ba- 
sic EFIT algorithm to include the fully non-linear 
toroidal flow and anisotropy contributions to the 2- 
D plasma equilibrium problem. The physical model 
is based on the guiding-centre plasma (GCP) for- 



malisr 



12j as derived by Dobrott and Greene |13| 
for a two-temperature anisotropic plasma model |2j . 
This new algorithm has been demonstrated in EFIT 
TENSOR, a code created by modifying an existing 
EFIT implementation used currently for the MAST 
tokamak (known as 'EFIT++' [3|). We will present 
analytical and MAST test case equilibria produced 
by EFIT TENSOR and demonstrate correct numer- 
ical convergence to force balance. 



2. Basic equations and assumptions 

We are concerned with the system of macroscopic 
equilibrium equations, based on the guiding-centre 
plasma and the MHD Ohm's law, given in natural 
units as [l2, 13 1 
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The single fluid pressure dyad P and its compo- 
nents pihPX) the single fluid velocity u and sin- 
gle fluid mass density p have textbook definitions 
(e.g.: [H|) in terms of each individual fluid equa- 
tion for each species present in the plasma, which 



in turn are moments of the guiding centre plasma 
model for each species. Divergence-less quantities 
(Eqs. © and in a 2-d axisymmetric cylindrical 
system (R, 0, Z) permit the covariant representa- 
tions 



B = Vt/> x V0 + RB<f,V(t> 
pu = VtpM x V</> + Rpu^Vcp 



(8) 
(9) 



in terms of poloidal stream functions tp and ipM- 
Combining Eqs. © and © with Eqs. © and © 
gives the well-known 'frozen-in' condition for mag- 
netic field lines in axisymmetric cylindrical geome- 
try 



B + R 2 VL{ip)Vcj) 



(10) 



which specifies that, to conserve magnetic flux 
through a given fluid element, all macroscopic ideal 
flow must occur parallel to the magnetic field or in 
a symmetry direction. Under these assumptions, 
flow in a poloidal direction can only manifest as 
the projection of flow in the parallel direction and 
toroidal angular velocity Q is constant on poloidal 
flux surfaces. In this study, we neglect all poloidal 
flows setting ip' M {ip) = 0. 

To close the system of equations, an energy equa- 
tion is required. Ignoring resistivity or other dissi- 
pation, the work done against the pressure equals 
the change in plasma energy U for a reversible pro- 
cess 

*j> ^ P\\ dp AB dB dU (II) 
p dt dt dt 

(noting that under most circumstances in this work, 
we will keep all thermodynamic quantities in units 
of energy density per unit mass). In the guiding 
centre model, relationships between the moments 
of the distribution function are obtained from a ki- 
netic equation (l2j |. These relationships depend on 
the form of the distribution function rather than the 
macroscopic variables alone, which is inconvenient 
for a fluid description. However, a thermodynamic 
statement can be reconciled @, [l5[ with the kinetic 
model and give the same results for certain simple 
choices of the functional form of the macroscopic 
variables. For example, an isotropic Maxwellian 
distribution, p = pT(tp), gives identical results in 
both models. In this study, our GCP compatible 
choice is to assume a two-temperature Maxwellian 
distribution of the form 



P\\(p,B,tp) = pT\\(ip) 
_{p,B,il>) = pT ± (B,1>) 



(12) 
(13) 
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Grad-Shafranov equation with toroidal 
flow and anisotropy 



Many examples of the Grad-Shafranov equation 
exist for static 



flowing [It 




anisotropic 
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21( equilibria. An equation in- 



coporating both flow and anisotropy @, |l5[ is used 
in this study and is outlined in this section. Vari- 
ous reductions of this general case to other known 
forms are given in the appendix. 

The basic scalar equations for force balance are 
obtained from components of Eq. (fTJ) . The toroidal 
component of force balance yields a new flux func- 
tion F(ip) for the toroidal magnetic field 

Fty) = (1 - A)RBi (14) 

the parallel component of force balance gives a 
Bernoulli equation and new flux function H(ip) 

H{^)^W{p,B^)-^R 2 ^f (15) 



W(p,B,^) = U + 



P 



(16) 



written in terms of a Legedre transform from energy 
U to enthalpy W in Eq. ([16)1. Finally, we recover a 
modified Grad-Shafranov equation from the ip com- 
ponent of force balance 
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subject to the integrability conditions constrained 
by energy conservation (Eq. (|lip) 
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Substituting the two-temperature Maxwellian GCP 
plasma expressions (Eqs. (|12[) and (TT31) ') into the in- 
tegrability conditions (Eqs. (|T8)) and (fH?)) ) gives a 
choice for enthalpy Wqcp(p, B,ip) which is consis- 
tent with with the GCP theory 
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for arbitrary choice of po and H gaugc . The gauge 
transformation is possible because of the freedom to 
choose (dW/dip) B and is physically related to set- 
ting an arbitrary reference energy in the Bernoulli 
equation (Eq. (|15p ). With the above assumption 
about parallel heat transport, the system of equa- 
tions is closed, and Eq. (|T7)) is now a second-order 
partial differential equation fully specified by the 
five flux functions: 



(22) 



with definitions coming from Eq. (fl2|) . Eq. (ITSj) . 
Eq. 00]), Eq. O and Eq. (EH). The tokamak 
equilibrium problem has been reduced to solving 
Eq. (fTTJ) by fitting the flux functions (Eq. ([2"2|)) to 
experimental data for appropriate boundary condi- 
tions. 

4. EFIT TENSOR 

4.1. EFIT 

Codes based on the original EFIT reconstruct 
the toroidal current profile from experimental mea- 
surements or given values, assuming an MHD force 
balance parameterisation for the current [l[. The 
Grad-Shafranov equation for isotropic and static 
cases is given by 



A*V> = i? 2 V 



I R 2 



= -R 2 p'{ip) - FF'(ip) (23) 



which is the limiting case of the more general sys- 
tem (Eq. (TIT])) for Cl(ip) = and Q(ip) = (al- 
though toroidal flow can be approximately incorpo- 
rated into 'pressure' in some codes, see appendix). 
The covariant representation for the field (Eq. ©) 
and Ampere's law (Eq. (|SJ)) give the identity 



A> 



-RJ<b 



(24) 



which is true for any 2-D axisymmetric equilibrium. 
This implies a parameterisation for the toroidal cur- 
rent 

J^R, V) = -Rp'W - FF'(i;)/R (25) 

when isotropic and static force balance is assumed. 
The two flux functions are decomposed into a linear 
combination of basis functions with constant coef- 
ficients 
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where the b(ip) represent the basis functions which 
are a priori assumed. The EFIT++ code currently 
supports polynomial, tension spline or Chebyshev 
polynomial representations. Eq. (j24|) is solved by 
alternately fitting to experimental or modelled 
constraints then performing a fast Buneman inver- 
sion (H of 



A*(^™ +1 )) = -RJ 4> (R,^) 



(28) 



at each nth iteration. The resulting ip(R, Z) 
solution is completely consistent with the input 
Jcf,(R,Z), but a recalculation of the force balance 
criterion Eq. ((251) changes J<p(R, Z) from the pre- 
vious iteration. Thus, a self-consistent Picard iter- 
ation occurs between if; and J^R.ip) in Eq. (|2"8")) 
until the change in ip(R, Z) is below an arbitrary 
threshold. The free flux functions are fitted to the 
experimental and given constraints with error by 
minimising 
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where the N are the number of data points, the a 
are the uncertainty, the C are the calculated values 
and M and P are the measured and given values 
respectively. In addition to the plasma currents in a 
tokamak, other specified, induced or unknown cur- 
rents are present. When including these known, un- 
known and plasma currents, the calculated poloidal 
magnetic field at a given magnetic probe location 
is a superposition of current contributions through 
a linear combination of the current coefficients 
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A;— free 

G(r h r)J^(R,^)dRdZ 

/ plasma 

where G is the response function for the probe. 
Other constraints are also expressed as a linear 
combination of the flux functions. Thus, the free 
coefficients and unknown currents constitute a lin- 
ear least-squares problem expressible as 



X 



= \\Au-n 



(31) 



where u contains the unknown coefficients and cur- 
rents, and k contains the constraints. This stan- 
dard problem is solved in EFIT using singular value 
decomposition (SVD). 



4.2. EFIT TENSOR system of equations 

Here, we explicitly present the EFIT TENSOR 
system of equations in S.I. units 

A> = -V0RJ4, 
(32) 



+RpH'(tp) - R 
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(34) 
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(36) 

A _ P\\~P± 

(37) 

We have re-arranged Eq. (|17p into an expression in 

taking advantage of the very general Eq. 
We have also re-arranged the Bernoulli relation 
(Eq. (1151) ) into an explicit form for the mass density 
p. At fixed p, the current is almost a linear com- 
bination of the flux functions TLD.D,' , FF' and a 
non-linear function of O. This system of equations 
is subject to the following identities for the five free 
flux functions 



p\\ m 
p k 



F(if>) = RB (j> \l- 
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P \PoP±J 2 

ew = -pb (-- — 

m \p\\ p± 



(38) 
(39) 
(40) 
(41) 
(42) 



These identities are straightforward functions of the 
kinetic moments p, v<p,p±,p\\ and field at any (i?, Z) 
location. The MHD particle mass m and Boltz- 
mann constant fc are only included to give Tji di- 
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mensions of temperature, but a more useful quan- 
tity is obtained when k = 1, m — 1 and Xji has the 
interpretation of the ratio of parallel pressure pn to 
mass density p. 

4-3. Numerical scheme 

The EFIT TENSOR code is a significant alter- 
ation to EFIT with a completely different set of 
physical assumptions, equations and free functions. 
However, many of the methods and constraints were 
adaptable to the more general case, and here we de- 
scribe those adaptations. 

The most important difference characteristic of 
the more general system is that (Eq. (J33J)) can- 
not in general be expressed as J<f, = J<j,(R,ip), nor 
can it be completely expressed as a linear combi- 
nation of free flux functions. The problematic con- 
tributions are due to A, p and (dW/dip) p B - We 
parametrize the linear current terms and also 0(V>) 
in terms of the basis functions b(ip) 



i=i 

i=l 

i=i 



(43) 
(44) 
(45) 
(46) 
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which are included in the iteration scheme 

j(n+l) _ R k (n) d T (n+1) 



1 



p R(l - A(")) dtp 
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(n+l) 



+#>) A_ff(»+i) _ RX^+^AtR. z)(») (48) 
ay 

where we have introduced a non-linear plasma cur- 
rent A(i?, Z) and weighting A corresponding to the 
bracketed term in Eq. (|3"3"]l. 



A(i?, Z) 



/cW\ 1 /A 
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This contribution is treated as a new basis function 
computed with bi-cubic spline derivatives with an 
associated nuisance coefficient A which is equal to 
unity for a converged solution. The reason in incor- 
porating this non-linear term in this way instead of 
defining a fixed current was to allow greater flexi- 
bility for intermediate iterations at the cost of an 
additional degree of freedom. Forcing A = 1 would 
be equivalent to specifying a known current based 
on the previous iteration. Contributions to the cur- 
rent from p and A in Eq. (|48l) are calculated from 
the previous iteration using Eqs. (|34p and (|3T| . 

It is clear from Eq. (|4"8")l that, in addition to the 
iteration of ip in Eq. (|25[) , the current which satis- 
fies force balance must also Picard iterate p, B and 
A contributions. Of particular concern is the expo- 
nential dependance of the density p on the flux func- 
tions H, Tii , and fi, which could conceivably pro- 
duce large variation in the inferred flux functions 
with each iteration if not suitably constrained. 

Eq. (|4"8"|) is a linear expression in terms of ba- 
sis functions which can be inserted into standard 
EFIT current constraints such as Eq. (I3TJ1) . In ad- 
dition, the flux functions may be independently 
constrained using the identities in Eqs. ((38)) - (|42|) . 
These identity constraints may be weak or strong; 
weak constraints depend on the kinetic moments 
p, or and the intermediate calculation 

for B whereas strong constraints also provide B. 
An additional constraint should be included which 
specifies that a = 1. The weighting on this condi- 
tion may also be specified to aid with convergence. 

An important feature of EFIT is the ability to 
infer the flux functions p'(ip) and FF'(ip) from the 
current profile. This inversion is possible due to the 
different R dependence in Eq. (|25|) . However, this is 
not possible with the more complicated expression 
Eq. dHJ). The degeneracy ofT!M) and H'ty), with 
each term weighted equally in R and p means that 
the columns of the linear inversion are linearly de- 
pendant and impossible to distinguish from current 
measurements alone. Furthermore, the anisotropy 
flux function 0(V>) has no linear contribution to 
the current at all, but instead influences the cur- 
rent through p and A. It is clear that EFIT TEN- 
SOR is a 'forward code' when invoking flow and 
anisotropy physics with additional information re- 
quired than the usual current constraints to resolve 
the degeneracy. This is the price paid for the fast 
linear inversion of an EFIT algorithm. A non-linear 
least-squares fit that includes density would break 
this degeneracy for significant computational cost 
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and complexity. 

To date, EFIT++ supports a number of con- 
straints, and those that have been adapted to EFIT 
TENSOR are listed here for reference: vacuum 
toroidal field, flux loops, magnetic probes, plasma 
current, poloidal field coils, safety factor on axis 
qo, static and rotational pressure approximations, 
B components, diamagnetic flux, boundary, equal 
if) and Motional Stark Effect (MSE). 

4-4- Dirichlet fixed-boundary mode 

The magnetostatic relation Eq. that is solved 
by EFIT at each iteration, constitutes a second- 
order inhomogenous partial differential equation on 
a finite discrete rectilinear domain with homoge- 
nous boundary conditions specified at infinity. For 
a given current distribution on that domain, the 
solution is unique. EFIT is a so-called 'free bound- 
ary' code which iterates the current distribution 
on this domain to give the best fit to the (usually 
measured) constraints and external currents. The 
external currents are critical to the shape of the 
plasma boundary on this domain. 

When considering analytical solutions, such as 
the Soloviev solution [23|], or other fixed-boundary 
problems, the external currents which give rise to 
the boundary are unknown. In particular, for the 
Soloviev case, the currents extend to infinity and 
cannot be included directly on a finite grid. For the 
purposes of benchmarking against an analytical so- 
lution, a fixed-boundary mode was added to EFIT 
TENSOR. This was done using a 'capacitance ma- 
trix' method (see for example which we will 
briefly describe here. 

The 'capacitance matrix' method is a way of 
calculating what surface currents are required on 
an arbitrary irregular boundary to satisfy a set of 
boundary conditions for ip on that boundary. Sup- 
pose we wish to solve Eq. (fM|) on a subdomain 
uj bounded by duj with inhomogenous Dirichlet 
boundary conditions ip(s) = f(s), Jj>(s) =0,sedw 
for some arbitrary /. If we only care about the solu- 
tion on the subdomain uj, we may solve the problem 
using the Green's function from the infinite domain 
fl using a sequence of steps: 

First, solve the inhomogenous equation for i\)\ on 
n with boundary conditions at infinity 

A*V>i = -RJ<p (50) 

and obtain a new function ipt> by measuring ipi on 



the boundary of the subdomain dw 

^6(s) = tpi(s),s e dw (51) 

Next, define a different homogenous problem for ip2 

A*V> 2 - (52) 

subject to the inhomogenous boundary condition 
using the function ipb measured before "02 (s) = 
f(s) —ipb(s). The function ip2(s) may be expressed 
completely in terms of unknown surface currents 
Jtf,{s) and the Green's function G(s' , s) on fi 

ip2 (s) = / J^is^Gis^^ds' 

= f(s) - Ms) (53) 

For the discrete problem, the solution for the vec- 
tor J^) may be expressed as the inversion of the 
response matrix G containing only the points on 
the boundary 

G-\f-M^J<t> (54) 

then the sum tp = ipx + ^2 uniquely satisfies the in- 
homogenous Dirichlet boundary conditions ?/>(s) ~ 
/ (s), Jj>(s) = 0, s 6 duj on w and the fixed boundary 
problem is solved. 

In the current implementation of EFIT TEN- 
SOR, the user specifies which grid points consti- 
tute the boundary and what the value of ip is on 
each grid point. Since the continuous boundary 
will not, in general, cross the grid points, the val- 
ues at the grid points may vary accordingly. The 
last closed flux surface search is also disabled and 
specified as the supplied boundary. Fig. [T] is an ex- 
ample fixed-boundary solution using this method. 
The ip contours inside the limiter region correspond 
to the Soloviev solution, and the contours outside 
the limiter are a by-product of the method and of 
no interest. 



5. Tests 

5.1. Extension of Soloviev solution 

The Soloviev solution is an analytical solution to 
the Grad-Shafranov equation when p(ip) an( i ^(VO 
in Eq. ((23]) are linear functions of ip. Adopting a 



G 



the Grad-Shafranov equation as 



N 




Figure 1: A Soloviev solution as reconstructed with EFIT 
TENSOR using a fixed-boundary mode. 



useful form used by [25|, the solution can be ex- 
pressed as 



ie(l-. 



I- -A [1 + era (2 + ear)] (J)' 



(55) 



a 



tp (56) 



R = ax + R (57) 

Z = ay (58) 

ps$)=l/[l-i>] (59) 

Fl$)=F a [l-i,} + RlBl (60) 



where {e, cr, r} control the shape of the tp solution 
and are known as the inverse aspect ratio, elliptic- 
ity and triangularity respectively. The field scal- 
ing is controlled by a and Bq. The linear profiles 
{ps, Fg} correspond to the static and isotropic flux 
functions and their gradients {p',F 12 } depend on 
the above parameters. To make the simplest pos- 
sible extension to flow and anisotropy, we re-write 



I R 2 



1 



(1-A) 2 
( dp 



p x =p x (R,B,i)) (61) 
Vv, log (1-A) = 
F{iP)F>{^) 



(1-A) \ d^ /B>R 
'dp± 
dR 



R 2 (l-A) 2 (62) 
= pRQ,{ipf (63) 



4>,B 

( ^Ei 
^ dB 



= -AB (64) 



4>,R 



If we wish to maintain the same tp geometry 
for our flow and anisotropy solution as for the 
Soloviev solution, then the partial derivatives of p± 
and F with respect to ip m Eq. (1621) must be the 
same, whilst also satisfying the additional condi- 
tions (Eqs. (|63p and (|64|) ). By inspection, this is 
achieved with the profiles 



p ± (R, B, 0) = l Pa n 2 R 2 - ^B 2 + aoPsW (65) 



A 



P\\{R,B,^) = - Po n z R 2 + -fB 2 + a psW (66) 

F 2 (tJj)=<j 2 F 2 ^) (67) 

where flo,ao = 1 — Ao,po a re constants with re- 
spect to R,B,ip. On substitution of these expres- 
sions into Eq. (|6"2")l . the logarithmic second term 
disappears and we re-obtain the ordinary Grad- 
Shafranov equation and the usual Soloviev solution. 
It is interesting to note that this implies that flow 
and/or anisotropy shear are the required properties 
to alter the current profile and magnetic topology 
of the analytical solution. 

There are advantages and disadvantages to this 
solution; although this solution exhibits the impor- 
tant de-coupling of magnetic surfaces and pressure 
surfaces, it does not do so with any respect for 
transport physics or thermodynamic assumptions. 
More specifically, one cannot write the pressures as 
P\\(p,B,ip) = pT||(-0) and p±(p,B,i/j) = pT±(B,tj}) 
as assumed in EFIT TENSOR as discussed earlier. 
For example, dividing Eq. (|65p by density po will 
give T± = T±(R, B, iji). This is interesting because 
it implies that the inclusion of flow and anisotropy 
in equilibrium introduces energy transport into the 
force balance inference problem. A related effect 
has been noted previously |2j when considering the 
direction of density shift for different enthalpy as- 
sumptions. It is, therefore, to be expected that 
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when using this analytical solution as a constraint, 
the physical assumptions of the parallel heat trans- 
port model will prevent the same solution being ob- 
tained. 

5.2. Comparison to analytical solution 

5.2.1. Force balance convergence 

A series of force balance numerical benchmarks 
were carried out on EFIT TENSOR by constrain- 
ing to analytical Soloviev solutions and testing force 
balance of the resulting solution. A 2-point finite- 
difference comparison between the left and right 
hand sides of the fluid force equation (Eq. (TT])) 
was used to measure force balance error. The lin- 
ear current profile of the Soloviev solution meant 
that this force balance test was found to be ac- 
curate to 1 x 10~ 14 when the analytical solution 
was tested with the same scripts. The parameters 
used for the analytical constraints are shown in Ta- 
ble [T] The ordinary tokamak parameters such as 
field strength, plasma current and geometry were 
chosen to resemble an ITER like configuration, and 
the kinetic properties of flow, mass density and 
anisotropy were chosen to maximize the changes 
to the analytical pressure from the ordinary case 
whilst still converging for the same numerical op- 
tions such as grid resolution and polynomial order. 
The results of the force balance benchmark are pre- 

Table 1: Extended Soloviev solution parameters 



flow cases, which is likely due to the spline deriva- 
tives used for the non-linear current calculation of 
V • (-^-V^). Better numerical evaluations of this 
term are perhaps possible, but have not been pur- 
sued in this work. 
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Figure 2: Comparison between plasma and field force bal- 
ance at the geometric axis for extended Soloviev solutions. 
The parameters used for these cases are shown in Table [TJ 
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sented in Fig. [5] The figure shows that the code 
accuracy scales equally well for static, flowing and 
anisotropic cases, with near identical accuracy for 
flow and static cases. The anisotropic cases were 
approximately 3 times worse than the static and 



5.2.2. Analytical solution reconstruction 

In this section, we compare the reconstructed so- 
lutions to the analytical solutions for the three cases 
of static isotropic, flowing and anisotropic plasma. 
All three analytical solutions were generated on a 
257x257 grid. Kinetic moments of p\\,p±, v$ and p 
were taken from the generated analytical solutions 
and used as input constraints for EFIT TENSOR. 
For the ordinary static case, the only free function 
was Xjj with n t = 1 the order of the polynomial. For 
the flowing case, H' was required for convergence 
and was set to nh = 2. For the anisotropic case, O' 
was required instead, also with order ng — 2. These 
were the minimal order polynomials found to con- 
verge total ij) to better than 1 x 10~ 15 , lx 10~ n and 
1 x 10~ 7 for the ordinary, flowing and anisotropic 
cases respectively. 

A comparison of reconstructed radial current and 
pressure profiles is given in Fig. [3] for each of the 
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three cases. The only constraint required to obtain 
the static case was total current, whereas the flow- 
ing and anisotropic cases had identity constraints 
for T h n,H and 6 using Eqs. ( IMf and (p]). No 
current or field profile data was used as an input in 
any case. The pressure profiles in Fig. [3] show very 
good agreement; in the isotropic case, this is due 
to the uniqueness of equilibrium for a given current 
scaling assuming a linear pressure form. 

In the flow case, pressure was constrained, so 
it is unremarkable that the pressure profiles agree, 
however what is remarkable is that the pressure is 
displaced whilst the magnetic surfaces remain un- 
changed, as expected for large flows (corresponding 
to sonic Mach number ~ 0.8 and Alven Mach num- 
ber ~ 0.2 — 0.8 across the midplane). Fig.@]is a plot 
of the EFIT TENSOR solution clearly showing the 
pressure displacement, whilst correctly satisfying 
force balance in Fig. [2J As the flow increases with 
radius, the agreement between the reconstructed 
current and the analytical solution diverges as the 
heat flow assumption T = T(ip) becomes less rele- 
vant to the analytical case. 

The anisotropic case shows similar unremarkably 
good agreement for the constrained pressure, but a 
significant change is observed in the current profile. 
The disturbance is symmetric in poloidal flux ip and 
is a direct consequence of the unphysical B 2 depen- 
dence on Tii in the analytical solution. Both the 
analytic solution and the EFIT TENSOR solution 
satisfy force balance but for different parallel heat 
flow assumptions in Eq. (|I7[) . This demonstrates 
that even with very good agreement for pressure 
profiles, the inferred magnetic topology can be rad- 
ically different when considering anisotropy (even 
between anisotropy models), a result which quali- 
tatively agrees with previous findings on q profile 
in the presence of anisotropy 



Conversely, 



If 



it follows that any pressure inferred from a mea- 
sured magnetic topology must be done so after 
careful consideration of the plasma transport and 
anisotropy to be physically relevant. 

5.3. MAST TRANSP constraints 

An equilibrium reconstruction was performed for 
MAST discharge 18696 at 290 ms using standard 
MAST EFIT++ magnetic constraints; magnetic 
probe, flux loops, field and induced coil currents 
and total plasma current. Motional Stark Effect 
measurements were unavailable for this shot. In ad- 
dition to the usual constraints, identity constraints 
on pu , p j_ , I),/, , p were provided from TRANSP [26| to 



N 




Figure 4: The displacement of pressure contours from mag- 
netic surfaces in the reconstructed flowing Soloviev bench- 
mark. Pressure shown in red. 



constrain the five flux functions, thereby including 
full anisotropy and flow in the force balance. The 
conversion from equivalent documented TRANSP 
quantities was performed as (ignoring the CGS unit 
conversion) 

v 4> = Rx OMEGA (68) 
/) = NIxm, + a.m.u. x AIMP x NIMP (69) 
PU = PPLAS+2xUFASTPA (70) 
p x = PPLAS+UFASTPP (71) 

where NI, AIMP, NIMP, OMEGA, PPLAS, 
UFASTPA and UFASTPP are: total ion density, 
average impurity mass and density, toroidal angu- 
lar velocity, total thermal pressure, fast ion parallel 
and perpendicular energies. These TRANSP val- 
ues were mapped to the mid-plane as functions of 
major radius. The resulting best fit is shown in 
Fig. [3] The largest discrepancy between TRANSP 
and EFIT TENSOR in this example is evident in 
the centrifugal force balance. Even though density 
and rotation were prescribed, an EFIT TENSOR 
solution with a lower rotation was inferred. This 
is likely related to the density being prescribed as 
a flux function in TRANSP, and the consequent 
density symmetry in radial mapping. The 2-D 
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Figure 3: Reconstruction of analytical radial profiles using EFIT TENSOR and constraining to pressure. The difference in 
parallel heat transport assumptions of the analytical solution and the two temperature GCP model has resulted in a different 
current profile for the same pressure profile. 
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TRANSP field variables BR, BZ and BPHI were 
extracted and used in a radial force balance bench- 
mark similar to that used in previously in this paper 
(Fig. [6]). The force balance test in Fig. [6] shows 
that, for this example, TRANSP underestimates 
the plasma pressure contribution either side of the 
magnetic axis and that MHD equilibrium with flow 
and anisotropy is not satisfied away from the mag- 
netic axis. A discrepancy is not surprising since 
TRANSP only uses a rotational pressure approxi- 
mation and takes beam pressure to be the average 
of parallel and perpendicular components, but the 
magnitude of the error is interesting. The addi- 
tional structure in the TRANSP force profile is due 
to the flat spots in the assumed TRANSP pressure 
profile in Fig. [5j which is not replicated with the 
second order polynomial used in this EFIT TEN- 
SOR example. 

We have thus demonstrated an equilibrium re- 
construction of a real MAST discharge including 
flow and anisotropy corrections to full order us- 
ing measurements of radial profiles taken from 
TRANSP. We anticipate, in future work, to replace 
many of these constraints with experimental mea- 
surements of density, rotation and temperature. 

6. Conclusion 

A new force model for EFIT has been pre- 
sented which includes full order toroidal flow and 
anisotropy effects under the assumptions of ideal 
MHD for a single fluid guiding centre plasma. The 
parallel transport model assumed is that parallel 
temperature is a flux function. This model has 
been demonstrated and tested in the existing code 
EFIT++ and was found to produce identical results 
in the static and isotropic limits with a compara- 
ble computational cost. The model was also tested 
for force balance using a finite difference scheme 
and was benchmarked against an analytical solu- 
tion. It was found that significant differences in the 
inferred current are possible for the same kinetic 
plasma constraints which result from the choice of 
parallel heat transport model. This implies that for 
sufficiently anisotropic plasma, the parallel trans- 
port model can be compared against measured cur- 
rent profiles, providing a novel measure of heat flow 
from equilibrium constraints. 

In future work, we will investigate the effects of 
flow and anisotropy on MAST shots using EFIT 
TENSOR. We anticipate that rotation and den- 
sity constraints will be measured from experiment, 
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Figure 5: Reconstructed radial profiles for MAST18696 at 
290ms. 
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Figure 6: Finite difference radial force balance check for profiles produced by TRANSP and EFIT TENSOR, with the inclusion 
of full order flow and anisotropy, for MAST discharge 18696 @ 290ms. 
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and parallel and perpendicular pressures will be 
calculated using an NBI particle model such as 
NUBEAM US or LOCUST \M- 



7. Appendix 

7.1. Reduction of modified Grad-Shafronov equa- 
tion to isotropic and/or static cases 

It is instructive to consider the reduction of the 
more general equilibrium problem to simpler cases, 
and the corresponding reduction in the number 
of free parameters. Indeed, treatments including 
poloidal flow have a free flux function in addition 
to the five presented in this work (but, it should 
be noted, for a different expression of W(p, B, ip) 

am). 

Setting A = in Eq. (|T7]) gives the isotropic force 
balance problem for toroidal flow 

\R 2 
-pT'^j) - pH'ty) 



tion 



\dip J 



R 2 



//(D^inln^J-^) 2 (73) 



W{p,i } )=T{i>)\n\j-\ (74) 
{TW),H(i;),SlW),F(i>)}(75) 



which is the form found by many authors [18l . ll9l . l2g . 
l2ll ] . A consequence of Eq. (|75|) is another widely 
known and useful form for this system 



V ■ 



(76) 



V</>\1 fdp\ F(tp)F'(ip) 

p(R,i>) = p{R^)T{i>) (78) 
{p(R,iP),n^),F(^)} (79) 



If we further assume zero flow by setting Cl(ip) = 
in Eq. (fT7]) then, on re-examination of Eq. (|73p , we 
identify no more explicit radial dependence and a 
redundant separation of temperature and density. 
Thus, we arrive at the basic Grad-Shafranov equa- 



{ R 2 



= -p'W 

H{ip) = T(tp) In 



R 2 



(80) 



(81) 
(82) 



Po 

P (i>) = p{^m) 

{p(ip),F(ip)} (83) 
Conversely, the purely anisotropic system becomes 

v.[ (1 - a ,(f); 

-pT^-pH'^) 
fdW\ F(ip)F'(ip) 



\ dip J 



B,p 



ff(^) = T||(V)ln 
W{p,B^) = T^)]n 
T±{B,f/>) 



R?{1 - A) 

P W) 

p T ± (B,i>) 

P T0) 
PoT ± (B^) 

BT0) 



(84) 
(85) 
(86) 
(87) 



|S-7]|(V>)eW0| 

{T^),H^),F{^)^)} 
and analogous to the pure flow case, we perform the 
full if) derivative of the Bernoulli equation and using 
the integrability relations (Eqs. (|18l) and (1191) ) we 
obtain a similar set of equations 

' d P\\ 
dip 



V • 



(1 



^4 



F{^)F'{iP) 
i? 2 (l-A) 

P\\ {B, ip) = p(B, ijj)T\\ (ip) 

{pii(s,V),ew,F(v)} 

which reduces immediately to the form 
GradfI3. 



(89) 

(90) 

(91) 
(92) 
riven by 



7.2. Rotational pressure approximation 

Here we relate our set of free functions (Eq. (|22p) 
to the 'rotational pressure' approximation used in 
some existing codes. Working in S.I. units, the def- 
inition of rotational pressure is through an expan- 
sion in major radius 



R 2 



(93) 



(94) 
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where Rq is some origin (ideally, the magnetic axis). 
Casting the isotropic pressure in terms of x and 
Taylor expanding around x = gives the required 
definitions 

k B 
P\\(x,B,ip) = pq t— — — , 
m \B - Q(ip)T\\(ip)\ 

' mH(ip) 



xTu(ip) exp 



-Daxis(V') = P0- 



x exp 
B 



feT||(V) 
mRl(x + l)n(ip) 2 

' jnH(ip) 



B-e(V0W) 



exp 



x exp 
k 



feT||(V) 

2fcT||(V>) 
T||(V>)L> axis (V>) (97) 



[5: 

(95) P 
(96) 



1 



P rot (V>) = ^R 2 n(^) 2 D^ is {^) (98) 

which is clearly only strictly appropriate for (i/j) = 
0, i.e.: when the system is isotropic. 
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